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ABSTRACT 

Various astrophysical studies have motivated the investigation of the transport of high energy particles in magnetic turbulence, either 
in the source or en route to the observation sites. For strong turbulence and large rigidity, the pitch-angle scattering rate is governed 
by a simple law involving a mean free path that increases proportionally to the square of the particle energy. In this paper, we show 
that perpendicular diffusion deviates from this behavior in the presence of a mean field. We propose an exact theoretical derivation of 
the diffusion coefficients and show that a mean field significantly changes the transverse diffusion even in the presence of a stronger 
turbulent field. In particular, the transverse diffusion coefficient is shown to reach a finite value at large rigidity instead of increasing 
proportionally to the square of the particle energy. Our theoretical derivation is corroborated by a dedicated Monte Carlo simulation. 
We briefly discuss several possible applications in astrophysics. 
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1. Introduction 

The scattering and the spatial diffusion of high energy particles 
off magnetic turbulence play a crucial role in many fields of as- 
trophysics. For instance they are key ingredients of Fermi ac- 
celeration processes because they directly control the efficiency 
and the rate of particle acceleration. They determine the prop- 
erties of the confinement of astrophysical objects from jets to 
galaxies and clusters of galaxies, and governs the transport of 
the particles through interplanetary, interstellar, or intergalactic 
space. Diffusion has long been described by a quasi-linear the- 
ory approach (Jokipii 1966, 1973), which allows us to calculate 
the diffusion coefficients when the turbulent field is significantly 
weaker than the background field. However, in many circum- 
stances the level of turbulence turns out to be large so that this 
standard picture requires extension. Several studies have exam- 
ined the transport properties in strong turbulence by means of 
numerical simulations, e.g. Giacalone & Jokipii (1999), Casse 
et al. (2002), Candia & Roulet (2004), and Fatuzzo et al. (2010). 
Most of these investigations have focused on the situation in 
which a large-scale turbulence cascades toward small dissipative 
scales - as in the Kolmogorov scheme - and in which particles 
interact with gyroresonant modes of the turbulence spectrum. 
From the point of view of the particle, the turbulence therefore 
occurs on large scales, as the coherence length of the magnetic 
field corresponds roughly to the maximal scale of the turbulent 
spectrum. 

However, in a variety of physical situations, the Larmor ra- 
dius of the particle can exceed the coherence scale of the tur- 
bulence. The transport of particles downstream of a relativistic 
shock wave provides a clear example of this situation. The mean 
field is there mostly transverse to the flow because of Lorentz 
tranform effects and shock compression, and the turbulence that 
is excited in the shock precursor is generated on microscopic 
plasma skin depth scales. In this case, perpendicular diffusion 
at high (possibly very high) rigidity controls the transport of the 



particles back and forth from the shock. More generally, the high 
rigidity regime likely plays an important role in the deconfine- 
ment process of particles of high energy, when their Larmor ra- 
dius exceeds the size of the astrophysical system. However, this 
high rigidity regime has received little attention so far, except for 
the pioneering study of Shalchi & Dosch (2009). The pitch angle 
scattering rate is known to increase in proportion to the square 
of the particle energy in this limit, but the behavior of the trans- 
verse diffusion coefficient, which is crucial in the above contexts 
deserves a careful analysis. This analysis is the objective of the 
present paper. It will be found in particular that even a weak 
mean field, as measured relatively to the turbulent component, 
can affect the scaling of the perpendicular diffusion coefficient. 

The present paper describes both a theoretical and a numer- 
ical study of diffusion at high rigidity. The theoretical aspects 
are discussed in Section [2] while the numerical simulations are 
presented in Section [3] Finally in Section |4] we summarize our 
results and discuss some applications. 

2. Transport of high rigidity particles with a mean 
field 

2.1. Notations and summary of previous results 

The transport of particles in magnetostatic turbulence is charac- 
terized by the reduced rigidity p, the level of turbulence 77, and 
the power spectrum of magnetic fluctuations in three dimensions 
(hereafter 3D) $3d(k). These quantities are defined as 

where ?l denotes the Larmor radius of the particle in the total 
(mean B and turbulent 8B) field B where B 2 = B 2 + SB 2 , 
e the energy of the particle, and ( c the coherence length of the 
fluctuations. 
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The turbulence level 77 is defined as 



(SB 2 ) 



(6B 2 ) + B 2 



where 77 — > corresponds to weak turbulence and 77 - 
sponds to pure turbulence with no mean field. 

The correlation function C(r) of the random field 



(2) 



1 corre- 



C(r) = 



(6B(x + r)8B{x)) 
{SB 2 ) 



(3) 



can be written in terms of the one-dimensional power spectrum 
S(k) oc£ 2 S 3d (k) 



C(r) 



/ dkS(k)sm(kr)/(kr) 
J dkS(k) 



(4) 



Casse et al. (2002) defined the coherence length as the scale at 
which C(r) is maximum; if the power spectrum takes the form of 
a broad-band truncated power-law S (k) oc (k/k m i a yP for k m [ n < 
k < k max and zero otherwise, one finds for the coherence length 
{ c - 0.77£~| . Alternatively, one can define the coherence length 
as we do here, to be 



t c = f 
Jo 



drC(r) , 



where one then derives in a straightforward way 



-Hp 

2 T] Jo 



dkk' l S(k) , 



(5) 



(6) 



and the presence of 1 /77 results from our choice of normalization 
for the power spectrum 



f 

Jo 



dkS(k) = T]. 



(7) 



where in practice the spectrum is bounded between k m i n and 
k m . dx . Both definitions for £ c coincide to within a factor close 
to unity. As a function of the spectrum index /3, the coherence 
length is close to either £ m j n on larges scale for /3 > 1, or to 
on small scales for/? < 1. 

The scattering frequency v s is defined as the reciprocal of 
the decorrelation time of the pitch angle of the particle, the lat- 
ter being defined relative to the direction of the mean field. As 
discussed in Casse et al. (2002), the scattering frequency can be 
written 



Jtt>l k ~ lsik)dk 



j S(k)dk 



(8) 



an expression that extends to the strong turbulence regime the 
results of the quasi-linear theory. This leads to the scalings 



pp- 2 (p«l) 



2 c 

(p>>1) - 



(9) 



The Bohm scaling holds only in the very special case where /3 = 
1. In addition to these quantities, the notion of correlation time 
also plays an important role because it measures the time beyond 
which a particle experiences a force that is decorrelated from the 
initial one, along the particle trajectory. It is then defined as 



Jo 



C(\Ax(T)\)dT , 



(10) 



where Ax(t) represents the displacement after a time r in the 
turbulence. In quasi-linear theory, only the unperturbed trajec- 
tory is inserted into this definition, although one can extend that 
definition with a diffusive trajectory as we later indicate. 

If a relativistic particle travels over a coherence length of the 
turbulent field without having displayed any wiggle, correspond- 
ing to the regime p » 1, then t c ~ i z jc. This correlation time 
is much shorter than the scattering time v^T 1 ~ rf x f?l z jc in this 
regime. The correlation time r c can be recovered from Eq. ( [10} 
by using the ballistic approximation Ax(r) ^ cr, which is ap- 
propriate in this regime p » 1, in which case Eq. Q leads to 
t c = (?r/2)(/3 - l)/3 _1 /t m | n /c ~ ejc. We note that in the spe- 
cial case where the power-law index of turbulence /3 = 1 (Bohm 
regime) Eqs. ^ and ^ lead to i c = (/i m in/4)log(/l max /A mill ), 
where /l m ; n and A max are the shortest and the longest wavelengths 
of turbulence. 

If a particle experiences a chaotic motion on a length-scale 
smaller than { c , corresponding to the regime p <K 1, then the esti- 
mate is more complicated to obtain but one finds that r c ~ pfitdc 
as follows. Since the correlation time remains shorter than the 
scattering time, Casse et al. (2002) proposed a heuristic estimate 
in which decorrelation arises out of the small-scale modes with 
wavenumber k > k m [ n p~ l , which give rise to gyroresonant in- 
teractions with the particle of rigidity p. The modes with wave- 
lengths longer than the Larmor radius (i.e. k < A; m i n p _1 ) construct 
the field line to which the particle is attached, hence do not cause 
decorrelation on timescales shorter than the scattering time. The 
above correlation time is indeed shorter than the scattering time 
and increases with p. The heuristic estimate for p < 1 is con- 
sistent with quasi-linear theory when 77 <K 1 and with numerical 
results in the strong turbulence regime (Casse et al. 2002) can 
then be written as 



dkk~ l S(k) , 



(11) 



k>k mm p- 



which bears some resemblance to the case discussed before for 
p » 1, except that p explicitly enters the sine function, since 
one must now follow the orbit of the particle around the field 
line, and the integral is limited to k > A: m i n p _I for the reasons 
given above. The calculation then implies that t c ~ pPtdc as 
announced. The particle trajectory undergoes decoherence be- 
fore traveling ( c because of the large number of wiggles in the 
random field. 

Thus, except for 77 ~ 1 and p ~ 1 for which the correlation 
time becomes comparable to the scattering time, a Markovian 
theory of the scattering process is appropriate, even if the tur- 
bulence is strong, stronger even than the mean field. This is an 
essential key for the present discussion. 

Independently of the rigidity, the parallel diffusion coef- 
ficient is always given by D\\ = c 2 /(3v s ), even in the strong 
regime of turbulence. As for the transverse diffusion coefficient, 
in the strong regime at low rigidities, it does not follow a law 
similar to the quasi-linear result but is proportional to D\\ (Casse 
et al. 2002) because of the magnetic field line wandering that 
transmits parallel diffusion in the transverse direction. Casse 
et al. (2002) found in particular that D ± = rj 2 - 3 D\\ at small 
rigidities, which rules out the conjecture of Bohm's diffusion. In 
the next section, we discuss the transverse diffusion in the large 
rigidity regime. 
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2.2. Transverse diffusion at large rigidity 

As mentioned previously, in the large rigidity regime p » 1, 
the correlation time is (much) shorter than the scattering time, 
hence we expect to derive the parallel and transverse diffusion 
coefficients using a Markovian description of the trajectory. In 
particular when p » 1, the velocity changes by 1/p only over a 
correlation time. This implies that significant changes in the ve- 
locity occur on timescales that are much longer than the correla- 
tion time. Therefore we can assimilate the effect of small-scale 
fluctuations to a fully decorrelated white noise on the relevant 
timescales. 

To calculate the particle transport in a random field, one has 
to use the solution of the differential equation that governs the 
evolution of the particle velocity v 

= [fi + «2«] • v . (12) 

The quantities £2o and 6&(t) are rotation operators developed as 
linear combinations of the generators of the Lie algebra of the 
rotation group, L\, L 2 , L3 



(13) 



In detail, £2o = Q.oB' Li/B {) , where B' Q denotes the r'-th compo- 
nent of Bo and Qo = c/ rL.o the Larmor pulsation defined with re- 
spect to the mean field. With this notation, £2o ■ v = Qov xBq/Bq. 
The operator 6£l(t) is decomposed in a similar way as the gen- 
erators of the rotation group, and 6Q. = c/r L , where r L is now 
measured relatively to 5B. 

To solve the equation of motion, one uses an auxiliary vari- 
able u that is defined as 
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where 



v(t) = R (t) • u(t), 
R Q (t) = exp(fQo) ■ 



We then define 

fi(o = hay 

one finds that u(t) obeys 

d 



££2(0 • R (.t) 



dt 



u = £1(0 u . 



This equation is solved as 
u(f) = T exp 



f £2(0) df' 

Jo 



•u(0). 



(14) 
(15) 

(16) 
(17) 

(18) 



Because the operator in the exponent is time dependent, to pre- 
serve the exponential character of the solution, a time-ordering 
operator T has to be introduced, as we now explain. 

We note that u(0) = v(0), thus the solution for v is given by 



v(0 = /? (0 Texp 



Jo 



£2(0) dt' 



v(0) 



(19) 



The regular part of the field generates the regular rotation matrix 
Ro(t), while the exponential accounts for the effect of the tur- 
bulent part. The time-ordering operator T maintains the chrono- 
logical order of the products in the non-commuting Afe) in the 
expansion of the exponential operator, i.e. 



T£2(0) • £2(0) = £2(0) • *2(0) if h>t 2 , 
= £2(0) ■ £2(0) if h > h , 



(20) 



and so on for higher order products. Alternatively, the time- 
ordered expansion can be written as a Dyson series 



T exp 



f £2(0) dt' 

Jo 



= 1 + 



V do... dt n £2(0) ■ ■ • £2(0 



). 

(21) 

We now use the following theorem that holds for a Gaussian 
stationary random process in the white noise limit. As discussed 
in detail in the Appendix, this is a direct generalization to any 
Lie algebra of a well-known result for a scalar random process, 
with no other restriction than the white noise assumption 



T exp 
Texp 



f £2(0) dt' 

Jo 



dt 2 £2(0) £2(0) 



(22) 



Various properties of the turbulent field can be considered i.e. 
that is either isotropic with no helicity, isotropic with helic- 
ity, or anisotropic with rotation invariance in the transverse di- 
rection. All these cases can be easily treated, although we fo- 
cus on two relevant cases: (A) 3D isotropic turbulence and (B) 
two-dimensional (hereafter 2D) isotropic turbulence in the plane 
transverse to Bo, with 6B ■ Bo = 0. 

We define the projection operators fi ± on the plane transverse 
to Bo and the projection operator ft\\ along Bo. We now define the 
correlation function of the random rotation matrices (5£2 



<(K2(0)cf£2(0)> = (m i (h)6Q. j (t 2 )}L i Lj 



(23) 



where 

(6Q. i (h)6Q. j (t 2 )} = 2r c 6(h - h) 



(24) 

The scalars (<5Qy ) and (<5£2^> characterize the relative strengths 
of the turbulence in the parallel (to Bo) and perpendicular direc- 
tions. In particular, for 3D isotropic turbulence, (6Q^_) = 2(5Q.^), 
in which case the above correlator becomes proportional to the 
identity. Then, using the properties of ft ± , 7t\\ and the L,-, one 
finds 



<(K2(0)(«2(0)> = -2T c 6(h-t 2 ) 



<(5Q 2 >1 - <(5Q||> 2 7r N 



,(25) 



where (SO. 2 ) = (5Q.\\) 2 + (^Q 2 ). We note that the above correla- 
tion function holds for ££2, which should not be confused with 

£2, the latter being the quantity of relevance for calculating the 
transport properties, as expressed in Eq. (|22|. However, 



<£2(0)£2(0)> = e"""°<5£2(O)e" ( '^" ) " (y£2(O)>e' 2 " , (26) 

and, because |ft||,e' no ] = ^j.,e f£1 °J = 0, the correlation function 
for £2 is the same as that for (5£2. 



Using Eq. ( 22 1, one then finds the solution for v: 



<v(f)> = /?o(0-exp 



d? 2 <£2(0)£2(0)> 



v(0) , (27) 
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where the average is taken over the possible realizations of the 
turbulent field. This leads to 



1 



Mt)> = *o(0-exp 
Using the properties of jt\\ and ft ± , this can be rewritten as 
<v(f)> = ^oW-{exp[-fT c (<<5Q 2 )-<5Q2))]7r|| 

*iW0) ■ 



v(0). 
(28) 



f exp -trMSO l )--(SQ t A ) 
Therefore, one derives the general results 

(V||(0)V||(0) = exp [-tr e (<^ 2 > - <«^>)] <V||(0) 2 ). 
In the transverse direction, 

<v x (0) • v x (0> = exp 



exp 



(29) 



(30) 



x T v x (0)-^o-v x (0) 
-tr c (((JO 2 ) - ^((5Q 2 > 
x cos(Q f) <v x (0) 2 ) . 



(31) 



The last equality follows from developing the exponential Rq = 
exp Ufio)> noting that £2o = £2()i-3 for Bo oriented along z, L3 " = 

(-1)"tt x , L 3 2 " +1 = (-1)"L 3 , V(0) • ?r x • v x (0) = <v x (0) 2 > , and 
T v x (0)-L 3 -v x (0) = 0. 

The parallel D\\ and perpendicular D ± diffusion coefficients 
are directly obtained from the correlation functions of the veloc- 
ity components after averaging over the initial velocities 



Dll = d/<vn(0)v]|(0>, 
Jo 

>+oo 

dx(v x (0)-v x (0>. 



D x = 



(32) 



Using Eqs. ( |30| > and pi) , this leads to 
On 



1^ 

3 v s 



It Vi 

D ± = -c 2 



where 



3" v 2 + Q 2 ' 

v s = r c (<c5Q 2 > - <<5Q 2 >) , 
v x = r c (<<5Q 2 > - i<5Qi) 



(33) 



(34) 



These expressions for D ± are formally similar to the results of 
the so-called classical diffusion theory, although they are ob- 
tained here under different physical assumptions; in particular, 
a strong turbulence situation is assumed. 

In case (A), for 3D isotropic turbulence, (5Q 2 ) = 2(dCl 2 ) = 

2<50 2 >, so that 



2 2 2 77 c 
v± = v s = -{SO. >t c = 

3 3 p L £ c 



(35) 



One may note that the expression for v s matches that derived 
from a random walk argument for pitch angle diffusion. We also 
note that the above calculation for v s may be applied to the 
regime p <k 1, as long as the correlation time is shorter than 
the scattering time. This is true in the case of v s = (2/3)77p^ -2 , 
which is the standard quasi-linear theory result. The result for 
the perpendicular coefficient cannot, of course, be extended to 
the regime p <K 1, as the above calculation does not account for 
field line wandering. 

In case (B), for 2D transverse isotropic turbulence, (SQ. 2 ) = 

0, (6Q.I) = (6C1 2 ), hence 



2v ± = (6n 2 )T Q 



(36) 



This demonstrates that the transverse diffusion coefficients fol- 
lows the scalings, which we express here for case (A), i.e. 
isotropic turbulence 



D x =s -ci c 



1 

" 2 

B 2 



ctcfflri (l«p«B/B ) 



B. 



(B/B « p) . 



(37) 







The transition between these two regimes takes place at p ~ 
B/Bo =i 77^ 2 (1 - riY 1 ! 2 , corresponding to v s ~ Q.q. At larger 
rigidities, the perpendicular diffusion coefficient remains con- 
stant, while the parallel diffusion coefficient continues to in- 
crease as p 2 . 

This result is supported by the numerical simulation that we 
now discuss. 

3. Numerical simulation of the transport with a 
mean field for high rigidities 

3.1. Numerical set up 

A Monte Carlo strategy is adopted to measure the diffusion co- 
efficients by integrating a large number of particle trajectories in 
given turbulent magnetic field configurations. Averages are then 
performed and statistical values of the diffusion coefficients de- 
duced. The numerical set up is presented hereafter: we first dis- 
cuss the construction of the magnetic field, then the integration 
of particle motion from Lorentz-Newton equation, and finally 
the estimates of the diffusion coefficients. 

The total magnetic field is expressed as B — Bq + 6B as 
before. The regular field is oriented along z, and SB is assumed 
to be isotropic in the three dimensions. An algorithm similar to 
Giacalone & Jokipii (1999) is used to construct the turbulent 
component of magnetic field SB by summing over plane wave 
modes {N mo d) with turbulent wavelengths extending from L m - m = 
1 = 27i/k m!lx to L max = 2n/lc m i n , the power spectrum following a 
truncated power law between £ m j n and k m . dx . In detail 



6B(x) 



2^ G n (k n )i;„ cos (k n -x+B„). 



(38) 



With Fourier modes of amplitude G„, and wave vectors k n = 
jt-ek isotropically distributed, the unitary vector is perpendic- 
ular to k„ in order to ensure that V • dB-Q, and B„ e [0, 27r] 
represents the random phase. The power spectrum is normal- 
ized by the turbulence parameter rj introduced earlier such that 
(6B 2 ) — BqT]/(1 - rf). For definiteness, the mode amplitudes are 
constructed according to a Kolmogorov cascade with logarith- 
mic spacing between wavenumbers: G„ oc k n . We note that 
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the details of the inertial range of the turbulence are not impor- 
tant because we are interested in the scattering properties at large 
rigidities, when the particle Larmor radius ?l is larger than all 
turbulent length-scales. For a detailed presentation of the numer- 
ical turbulent magnetic field construction, the reader is referred 
to Section 2.B of Casse et al. (2002) and Section 3 of Giacalone 
&Jokipii(1999). 

Several tests of the dynamic range of turbulence L max /L m in 
and the magnetic wave-modes N mo d were performed. The main 
difficulty is that the scattering timescales increase as a square 
of particle rigidity. For large rigidities, it is thus difficult to 
preserve the accuracy with time when achieving particle dif- 
fusion together with a realistic magnetic field model. To de- 
velop a simulation that operates over a few scattering times, one 
needs to achieve an integration time of at least 100pFl/c, as in 
our simulations. One must also strike a compromise with the 
number of plane wave modes to save computational time; val- 
ues of order 200 - 300 have emerged as a satisfactory compro- 
mise between accuracy and calculation time. To save computa- 
tional time, and because the small scales of the turbulent cascade 
are of little influence, the dynamic range has been shortened to 
^min/imax = 0.1. Tests performed with a larger dynamic range 
have provided similar results; the highest accuracy is obtained 
when modes are concentrated on the largest scale. It is explained 
physically by the high energy particles interacting only with the 
largest magnetic structures. 

Particle motion is solved using the Lorentz-Newton equation 
of motion that preserves its energy, hence its Lorentz factor y 



dv 
dr 



g 
ymc 



v x (B + SB) 



(39) 



At this point, we define the numerical rigidity p' = 2^r L /L max , 
which differs from the previous physical definition by a numeri- 
cal factor of order unity, as discussed earlier. The exact relation 
between i c and L max /2n depends on the dynamic range and the 
power-law index of turbulence. In the following, the conversion 
factor between both rigidities is derived using £ c 0.1L max , a 
good approximation for a Kolmogorov-type spectrum. 

The numerical integration of Eq. ( 39 1 is performed using a 
Bulirsch-Stoer schema (Press et al.1986). Once a large number 
of particle trajectories were calculated and stored, statistical av- 
erages instead being performed. Given the number of particles 
N p for each field realization and the number of field realizations 
^Vfieid, the diffusion tensor coefficient (z, j) is evaluated as 



Duff) = 



1 ^ (Xjjt) - x,(t ))(xj(t) - Xj (to)) n , k 



2N p N mi n £ 



t-t 



2At 



(40) 



The average is performed over different particle trajectories 
and different field realizations. For each value of p, we take 
A^fidd xJVp = 10 3 different trajectories with random initial ve- 
locity directions. The asymptotic value for t — > oo (plateau) 
is roughly constant and defines the actual diffusion regime. It 
gives the diffusion coefficient as Djj as t — > oo, precisely when 
t » VjT 1 . This method of coefficient estimation appears precise 
enough for an integration involving 10 3 particles. A complemen- 
tary technique consists of evaluating time correlations between 
velocities over particle trajectories. With 1000 particles in the 
transport regime studied here, this method is affected by numer- 
ical noise for the velocity correlation function, hence is not pre- 
sented. 



10000 



1000 



100 



0.1 



10 

p' = 2%!L I L 



Fig. 1. The diffusion coefficient variation is plotted in units of 
cL max /(27r) as a function of rigidity p' in pure turbulence (Bo = 
or T] = 1). The dashed line is drawn as a reference for a scaling 
Di S0 oc p' 2 . For p' > 1, D iso is indeed proportional to p 2 . 



Two different cases were investigated numerically: a pure 
turbulence situation (Bo - 0) an d a strong turbulent case with 
SB » Bo. Results are presented in the following sub-sections. 

3.2. Pure turbulence B = 

These simulations were performed to test the correctness and ac- 
curacy of the code. On theoretical grounds (see the appendix of 
Casse et al. 2002, Aloisio et al. 2004, Pelletier et al. 2009) and 
previous numerical works (Parizot 2004), we expect the diffu- 
sion coefficient to evolve as the square of energy (e.g. rigidity) 
when p' » 1 . 

Here we set rj = 1 and 6B isotropically distributed by con- 
struction, so that the three space directions are equivalent. The 
equivalence of the three directions was numerically verified in 
our simulations. The diffusion coefficient is evaluated as 



(Ax 2 ) + (Ay 2 ) + (Az 2 ) 



6At 



(41) 



Figure [T] shows numerical values calculated for p' going from 1 
to 100. The diffusion coefficient is plotted in units of cL m . dx /(2n) 
as a function of rigidity p'. A power law is observed for 1 < p' < 
100, as predicted by the theory, namely D, so oc p' 2 oc e 2 . One 
may be able to discern a slight deviation at p' close to 100. This 
purely numerical effect disappears when taking a larger num- 
ber of magnetic wave-modes on scales close to L max by defining 
^min/imax ~ F We retain however the current field configura- 
tion, taking this effect into account when interpreting the results. 



3.3. Weak mean field B < SB 

We now consider the case where a constant weak mean field Bq 
along the z direction is present. In this case, two different diffu- 
sion coefficients are defined D\\ = D zz and D ± - (D xx + D vy )/2. 
Overall, we explored five different levels of turbulence rj = 
{0.5,0.9,0.99,0.999,0.9999}, spanning five orders of magnitude 
in SB 2 /Bq. The rigidity p' ranges from 1 to 100 for each value 
of rj. At each calculation point {77, p], the coefficients are eval- 
uated by averaging over 10 3 particles (10 particles x 100 field 
realizations), as before. 
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Fig. 2. The parallel diffusion coefficient D\\ plotted in units of 
(cL m . dx /{2n)) as a function of p' for different degrees of turbu- 
lence 6B 2 /B 2 ) e [1, 9999]. Here D\\ °c p 2 , as in the case of purely 
isotropic turbulence without mean field. As long as 6B 2 /B 2 ) » 1, 
the strength of the turbulence does not influence the normaliza- 
tion of D\\. 



Fig. 3. The transverse diffusion coefficient D ± plotted in units of 
cL m . dx /{2n) as a function of p' for different degrees of turbulence 
5B 2 jB\ e [1,9999]. The diffusion coefficient saturates at p' ~ 
SB/Bq. Below this value, its behavior is similar to the parallel 
diffusion coefficient. Beyond p', its value becomes independent 
of particle rigidity. 



As shown in Fig. [2] the parallel diffusion coefficient retains 
the same dependence on rigidity as in pure turbulence, D\\ oc p' 2 . 
For r] > 0.5, the turbulence level has almost no influence on the 
value of D\\ - QScf^p' . Therefore as expected, the mean field, as 
long as it remains weak enough, has no influence on the diffusion 
of particles along its direction. 

The picture is different for the transverse coefficient when 
the particle rigidity becomes large. In Figure [3] the simulated 
transverse diffusion coefficient is plotted as a function of rigid- 
ity p' for different degrees of turbulence. In each case, its value 
saturates to a constant value when p' ~ 6B/B0. This value be- 
haves proportionally to the turbulence degree; in detail, D ± - 
0.13c(L m a X /27T)SB 2 / Bq, in excellent agreement with our theoret- 



ical prediction from Eq.(37i. Individual particle trajectories re- 
veal a weakly perturbed helical path when p' >> 1. Therefore, a 
strong small-scale turbulence acts as a collection of small-scale 
scattering centers, each producing a small deflection. 

According to the theory, D\\ is the limit of a function 
c 2 g||(f)/3 as t — > 00, precisely as t > t s , the function being 



1 - <r"»< 

fl(0 = — ■ 

In a similar way D ± is the limit of a function c 2 g ± (t)/3 as t 
in addition to when t > t s , the function being 



(42) 



n 2 + v 2 



1 -, 



cos (£2o?) sin (C2o?) 



(43) 



The numerical simulation reproduces these types of be- 
havior, although the transverse evolution departs slightly from 
the above formula before reaching the scattering time r s . 
Nevertheless, the agreement between the theory and the numeri- 
cal simulation holds during the linear growth at the beginning of 
the evolution and when the evolution approaches the asymptotic 
behavior. The numerical results confirm the theory we proposed 
in the previous section for the asymptotic regime. The scatter- 
ing time is clearly the time beyond which spatial diffusion takes 



place. We can also note that there is a sub-diffusion regime be- 
fore the settlement of the transverse diffusion regime. 

The anisotropy ratio D x /D\\ can be seen in Figure B] as a 
function of p'. When the turbulence level 77 is close to I and 
p' is not too large, the transport appears isotropic D x /D^ ^ 1. 
At higher rigidities, its behavior follows the law oc p'~ 2 for all 
turbulence levels, illustrating the saturation of the transverse co- 
efficient and in agreement with the theoretical prediction. 

3.4. Comparisons with previous results 

Transverse diffusion at high rigidity, as far as we know, has been 
poorly studied in the literature. However, we can compare our 
results with several previous numerical and theoretical studies 
with different limits. 



, , t .vs. D 

5B /B = 99, n = 0.99, p = 100 




t/t, 



Fig. 4. Transition toward parallel and perpendicular diffusion. 
Before reaching its asymptotic value for t > r s , the transverse 
diffusion rate decreases as in a sub-diffusive regime. 
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Fig. 5. Anisotropy ratio D x /D\\ as function of p' for different lev- 
els of turbulence SB 2 jB 2 Q e [1, 9999], as indicated by the various 
symbols. The dashed line provides a guide for a p'~ 2 scaling. 



The seminal study of Giacalone& Jokipii (1999) focused on 
the propagation of mildly relativistic particles (E = lMeV to 
1 GeV) in the interplanetary magnetic field (SB 2 ~ B^). Their 
simulations provided results for p < 1 and 77 < 0.5. However, 
they also performed several simulations in which the particle 
energy and the coherence length remained fixed, while the tur- 
bulence level was varied. In particular, they examined the case 
rhfi/^c = 10 for moderate values of SB 2 /B 2 ) (Fig. 6 of their pa- 
per) in which D ± /D\\ is plotted as a function of A||/rL,o (A\\ denot- 
ing the mean free path in the parallel direction). By inspecting 
their figure, one can see that they varied SB 2 /B 2 ] from 0.05 to 
30. As a result, they found a classical scattering theory scaling 
but no physical explanation was proposed. Strictly speaking, the 
classical theory is valid only for weak turbulence {SB 2 <K B^), 
which clearly does not apply to those simulations. The present 
theoretical framework provides a clear explanation of this result, 
which we confirmed with additional detailed numerical simula- 
tions. It is found, for instance, that particles with large rigidities 
do not interact directly with the magnetic field lines but experi- 
ence an overall magnetic topology dominated by the mean field 
with "infinite" coherence length. As a result, the particles exe- 
cute regular orbits around Z?o and undergo random deflections 
on the coherence length-scale. 

The simulations of Casse et al. (2002) investigated weak as 
well as strong turbulence regimes where SB 2 /B 2 Q e [0.1,99]. 
An FFT algorithm was used to construct the magnetic field. For 
p' > 1, these authors found evidence of anisotropic scattering 
D±/D\\ < 1 for all turbulence levels. However, only three simu- 
lations points were computed in the high rigidity range and the 
estimate of the power law slope was inaccurate. Nevertheless, a 
reasonable agreement is obtained when comparing values of Dy 
and D ± with the present results. 

Parizot (2004) presented simulations of particle propagation 
in pure isotropic turbulence. The results in the regime rt » £ c 
leads to a diffusion coefficient with a quadratic scaling, D oc E 2 , 
in agreement with our results from Sec |3.2| 

Shalchi & Dosch (2009) derived an analytical expression 
for the diffusion anisotropy ratio D ± /D\\ in the framework of a 
non-linear guiding centretheory. They assume an isotropic tur- 
bulence 6B with a mean field Bo- No assumption was made 
about either the level of turbulence or about particle energy, 



Fig. 6. Ratio D ± /D\\ as a function of p' for 77 = 0.99, compared to 
theorectical predictions and other numerical simulations. Filled 
diamonds: our simulation results. Star symbols: results from 
Casse et al. 2002. Solid curve: present theoretical prediction with 
best-fit from simulations (see Fig. [2]). Dashed curve: present 
theoretical prediction with analytical D\\ = c 2 /(3v s ). Dot-dashed 
curve: analytical prediction from Shalchi & Dosch (2009), their 
Eq.(15). 



hence their result should be valid for any particle rigidity and 
turbulent field strength. An expression of D ± /D\\ [Eq. (15) in 
their work] that depends on two parameters was obtained. The 
first parameter corresponds to the ratio of the mean free path 
(A\\) along the mean field direction to the coherence length { c 
of the turbulent field. The second parameter is the turbulence 
level SB 2 IB\ = 77/(1 - 77). Shalchi & Dosch (2009) thus find that 
the transport becomes highly anisotropic, meaning D x /D\\ <K 1 
when A\\/-£ c » 1 and/or SB 2 /B 2 ) is not too large (see Figs. 1 and 
2 of their work). Therefore, our present conclusions agree with 
theirs, at least at a qualitative level. A detailed comparison would 
require us to define A\\ as a function of p', which could be done 
by using our results of Dy for which p' = (4m]/30) l l 2 {A\\/£ c ) 1 / 2 . 
With this substitution, we can directly compare their predictions 
to our results. In Fig. [6] we plot the ratio of diffusion coefficients 
as a function of p' from our numerical simulations and compare 
these results to both the predictions of Shalchi & Dosch (2009) 
and the theoretical model developed in Sec. [2] Good agreement 
is found between the simulation results (diamond symbols) and 
our theory (solid curve and dashed curve); however, the predic- 
tions of Shalchi & Dosch (2009) disagree with the numerical 
results, increasingly so as the rigidity increases. In particular, 
their analysis predicts a scaling with a slope -2.4 instead of the 
value of -2 observed here. Repeating the same comparisons for 
each simulated value of SB 2 IB 2 , we were unable to find agree- 
ment between the predictions of Shalchi & Dosch (2009) and our 
simulations; the predicted values always lie below the numerical 
results, with a different power-law scaling, comprised between 
-2.5 for SB 2 IB 2 = 1 and -2.4 for SB 2 IB 2 = 10 4 . At this point, 
it could be argued that our definition of A\ as a function of p' is 
inaccurate. However, on physical grounds, the scaling A\\ oc e 2 
when Pl » £ c remains robust. Therefore, the discrepancy be- 
tween the power law scalings should not be affected by uncer- 
tainties in the numerical prefactors. We think that the "guiding 
center" assumption in their work is questionable. 
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4. Summary and some astrophysical applications 



converted into electromagnetic turbulence, i.e. 



4.1. Summary 

Our investigation of the diffusion process in small-scale turbu- 
lence with a mean field revealed that, despite its smallness, the 
mean field plays a role in transverse diffusion because the scat- 
tering frequency decreases like e~ 2 , whereas the Larmor fre- 
quency decreases like e _1 . Instead of finding a single diffusion 
coefficient that increases like e 2 , we found an anisotropic diffu- 
sion with a transverse coefficient that reaches a limit value at 
large rigidities. The theory we proposed is based on a single 
assumption, namely that the correlation time is much smaller 
than the scattering time, which is valid for both small and large 
rigidities. The only regime where the theory fails is for a rigidity 
close to 1 and a high turbulence level; however, the interpolation 
is obvious. The theory allows us to derive a correct pitch-angle 
scattering rate and a correct parallel diffusion coefficient for ev- 
ery rigidity. It provides a transverse diffusion coefficient similar 
to the classical scattering theory formula, despite the arbitrary 
level of turbulence, which is a correct result for large rigidity. At 
low rigidity, the present theory is incorrect because it does not 
take into account the effect of field line wandering described in 
Casse et al. 2002. 



4.2. Particle transport in relativistic shock environments 

One major application of the diffusion theory in small-scale tur- 
bulence is the transport of supra-thermal particles in the vicinity 
of a relativistic shock. By crossing the shock transition, electrons 
and protons reach more or less the same characteristic energy 
(e) ~ y^nipC 1 as revealed clearly by particle-in-cell simulations 
(e.g., Sironi & Spitkovsky 2011). There is a single plasma fre- 
quency Wp* ~ w p i, where o) p j is the ion plasma frequency in 
the upstream or unshocked plasma. This length-scale character- 
izes the typical length scale of the microturbulence excited in the 
shock precursor, as transmitted downstream of the shock transi- 
tion and viewed in the downstream rest frame. The generation of 
short scale intense micro-turbulence is possible only at low mag- 
netizations of the upstream plasma (Sironi & Spitkovsky 201 1), 
where the magnetization parameter <x is here defined as the flux 
of magnetic energy crossing the shock over the flux of matter 
energy, cr = B ( 2 sin 2 6b I '4np u c 2 (where 6b is the angle of the 
background magnetic field with the shock normal, and p u the 
unshocked plasma mass density). However, this same level of 
magnetization also permits the efficient acceleration of particles 
through a first-order Fermi process at the shock front (Lemoine 
& Pelletier 2010, 2011). For larger magnetizations - the exact 
level depending on the shock Lorentz factor, see the above ref- 
erences - the Fermi process cannot develop because of a lack 
of efficient scattering in the microturbulence (Lemoine et al. 
2006, Niemiec et al. 2006, Pelletier et al. 2009). In brief, the 
development of the Fermi process hinges on the development of 
micro-turbulence, which itself requires (in the absence of exter- 
nal sources of turbulence) a sufficiently low magnetization level. 
The situation in which particles are accelerated is by far the most 
interesting as it should produce directly observable signatures, in 
the form of radiation and possibly neutrinos. 

The transport properties of these accelerated particles is then 
directly governed by the parallel and perpendicular diffusion co- 
efficients in the limit of large rigidity, as discussed above. We 
assume that the microturbulence has a typical length-scale close 
to 6* - c/ojpt and that a fraction eg of shock dissipated energy is 



(SB 2 ) 
8tt 



= 2e B rshPuC 2 



(44) 



where the rigidity of shock accelerated particles of energy e is 
given by 

-1/2 & e 



(45) 



Current simulations indicate values of eg ~ 0.01 - 0.1, hence 
p > 1 and all the more so at high energy. 

In this regime, the perpendicular diffusion coefficient that we 
discussed in the previous section becomes particularly relevant, 
as the mean magnetic field is mostly perpendicular to the shock 
normal in the downstream frame, since the transverse compo- 
nents (relatively to the shock normal) are increased by 2 V2y sri , 
while the parallel component remains the same as in the up- 
stream frame. Therefore, perpendicular diffusion at high rigidity 
plays an essential role in the transport of particles in the down- 
stream flow of a relativistic shock. 

We consider the diffusive behavior of particles in the down- 
stream rest frame. In this frame the shock front appears to move 
away with velocity V sri ock - c/3. Achieving Fermi cycles re- 
quires the particle to return to the shock front. The return time is 
then measured by identifying shock front with the particle mean 
displacements 

|*ret = ^j2D ± t Kt . (46) 

Therefore f ret = 18Z) ± /c 2 and Fermi cycles are possible until 
f ret is neither large nor too short. While the first case is con- 
strained by confinement in the acceleration site, the second one 
is related to the diffusive approximation that is valid only when 
f re t ^ t s . Using the second limit to constrain diffusive returns, 
one obtain D ± /D\i > 1/6, equivalent to v s > V5Qo when D ± 
is replaced by its expression from Eq. (33 1. Fiducial values for 
a relativistic shock in the interstellar medium provide an energy 
limit Ey lm ~ 10 I9 eV. This limit is somewhat irrelevant because 
f ret » R-acc/c at this energy, where R^/c is the shock dynam- 
ics timescale. Hence, the returns appear to be efficient when the 
condition t s < f re1 <K R- dC c/c is satisfied. Further investigation 
would require us to solve a kinetic equation taking into account 
acceleration, scattering and energy losses processes. Diffusion 
coefficients obtained in this work may be relevant to providing 
more realistic results. Previous works assumed Bohm diffusion 
or isotropic pitch-angle scattering. 

Detailed discussion of the performance of the relativistic 
Fermi process is beyond of the scope of the present paper and 
is left to future work. 

In certain astrophysical settings, the transverse diffusion may 
play a key role in the transport of particles upstream of a rel- 
ativistic shock, most particularly if the shock propagates in a 
wind with a dominant toroidal field at large distances. These cir- 
cumstances can be encountered in particular when a gamma-ray 
burst explodes in the wind of the progenitor, or at the termination 
shock of a pulsar wind. 

4.3. High-energy cosmic rays 

The above result about transverse diffusion has a broader appli- 
cation than Fermi acceleration at shocks, as it governs the con- 
finement properties of any relativistic flow containing a small- 
scale turbulence, where "small" is measured relatively to the 
Larmor radius of the test particles propagating in this flow. This 
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concerns in particular the propagation of very high-energy cos- 
mic rays in our Galaxy. Assuming a coherence length of in- 
terstellar turbulence ( c ~ 10 - 100 pc, a mean field inten- 
sity of 3pG approximately and a turbulent field of the same 
order, the rigidity of particles of energy E is given by p 
2(£'/10 17 eV)(4/10pc)" 1 (B/5^G)" 1 , while the Larmor radius 
r L 20pc(£/10 17 eV)(B/5//Gr 1 . Assuming rj 0.5 in the 
Galaxy and using Eq. (33 I, the perpendicular mean free path is 
then of order A ± ~ 6 pc with these values of energy and magnetic 
field. This implies that the escape, or transport across the disk 
magnetic field of particles of energy > 10 17 eV is governed by 
the perpendicular diffusion in the high rigidity regime discussed 
above. Quite interestingly, this energy range presumably corre- 
sponds to the transition between the Galactic and extragalactic 
cosmic -ray components in the all-particle spectrum. 

Finally, one could mention another application of the present 
discussion, to the field of magnetic reconnection. There, trans- 
verse diffusion likely plays a role in the control of particle dif- 
fusion across the field lines with small-scale turbulence being 
associated with the dissipation of magnetic energy. The recon- 
nection rate depends on two fundamental parameters (Lyutikov 
& Uzdensky 2003): magnetization and the Lundquist number 
that involves diffusion across field lines. In general, one assumes 
Bohm diffusion for simplicity but the present work provides the 
grounds for a more accurate estimate. 

Appendix A: Average of the time ordered 
exponential 

We solve the differential equation u = fi • u in successive it- 
erations that leads to a Dyson series, the average of which is 
composed of products of the form 



rt pti r-hp- 

A 2p (t)= dh df 2 ... 

Jo Jo Jo 



(A.l) 

which can be compared to Eq. pi) . 

For a Gaussian process, each average of order 2p products 
can be divided into a sum of p products of second-order mo- 
ments, the sum containing (2p - 1)!! terms. We assume a sta- 
tionary random process such that the second order moment is an 
even function of the time difference. In the white noise limit, the 
"nested" and "crossed" averages vanish, only the "unconnected" 
averages remaining in the expansion. Nested terms contain prod- 
ucts of the form (X(tj)X(ti)){X(tj)X(t k )} with t t > tj > t k > t h 
while crossed terms are of the form (X(ti)X(tic))(X(tj)X(ti)) with 
£; > tj >t k > t[. These terms vanish as the various delta functions 
associated with the second order moments cancel each other as 
a result of the time ordering in the upper bounds of the integrals. 
Thus, only the unconnected average remains at each order 



dh df 2 ... dt 2p 
o Jo Jo 



(Slitt) ■ ft(f 2 )> • • ■ i&(t 2p -i) ■ W 2p )) 



(A.2) 



We introduce the short-hand notation (fi(fi) • fife)) = 
C(t\ — t 2 ) = 2T c 5{t\ - f 2 )CV Then one can calculate A 2p (t) by 
recursion, starting from the last double integral in the product 



A 2p (t) = 



(Tcty.p 
pi » 



(A.3) 



We consider now the integral of the second order moment 

K(t)= f d?! f df 2 (tlih) ■ h(t 2 )) = 2t c ?C„ . (A.4) 
Jo Jo 



Therefore 

Kit) = -Li»(o . 

hence summing all the terms of the series, 

p=+oo 



^ A 2p (t) = exp 



p=o 



(A.5) 



(A.6) 



Further details can be found in Frisch (1966) and Pelletier 
(1977). 
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